/*
 * @Description: 优化法估计矩阵1范数
 * @Autor: kaikaima
 * @Date: 2021-04-06 20:27:16
 * @LastEditors: kaikaima
 * @LastEditTime: 2021-04-29 19:24:19
 */

#include <iostream>
#include <cmath>
#include "L_equations.h"

int main(int argc, char* argv[]){
    for(int N=5;N<=20;N++){
        double** A;
        A=new double* [N+1];
        for(int i=0;i<=N;i++) A[i]=new double[N];

        for(int i=0;i<N;i++){
            A[N][i]=1.0/N;
            for(int j=0;j<N;j++)
            A[i][j]=1.0/(i+j+1.0);
        }

        double key1=0,key2=0;
        double *z,*w,v[N];
        L_Equations Hilbert(N);

        while(true){
            Hilbert.set(A);
            //Hilbert.transpose();
            w=Hilbert.Cholesky_solve_modified();

            for(int i=0;i<N;i++)
                if(w[i]>0) v[i]=1;
                else if(w[i]==0) v[i]=0;
                else v[i]=-1;
            
            Hilbert.set(A);
            Hilbert.set(v);
            z=Hilbert.Cholesky_solve_modified();
            
            //get ||z||∞ zTx
            int zn=0;
            double zx=0;
            for(int i=1;i<N;i++) if(fabs(z[i])>fabs(z[zn])) zn=i;
            for(int i=0;i<N;i++) zx+=(z[i]*A[N][i]);

            if(fabs(z[zn])<=zx){
                for(int i=0;i<N;i++)
                key2+=fabs(w[i]);
                for(int i=0;i<N;i++){
                    double tmp=0;
                    for(int j=0;j<N;j++)
                    tmp+=fabs(A[j][i]);
                    if(key1<tmp) key1=tmp;
                }
                break;
            }
            else{
                for(int i=0;i<N;i++)
                    if(i==zn) A[N][i]=1;
                    else A[N][i]=0;
            }

            delete[] w;
            delete[] z;
        }
        delete[] w;
        delete[] z;

        std::cout<<N<<"阶Hilbert矩阵的无穷范数条件数: "<<key1<<" * "<<key2<<" = "<<key1*key2<<std::endl;

        for(int i=0;i<=N;i++)
        delete[] A[i];
        delete[] A;
    }
    return 0;
}